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ABSTRACT 

A theoretical study is made of the shearing flow over a 
sinusoidal boundary representing the interface between the ocean 
and the atmosphere, with the main purpose being to investigate 
the pressure distribution at the interface and to calculate the 
energy transfer between the two media. 

The theory is developed on a model of turbulent flow making 
use of Prandtl's mixing length theory to represent the shear 
stress terms in the basic Navier-Stokes equations. A curvilinear 
coordinate system which follows the wave train is used in order 
to simplify the equations, and to allow for a linearized solution 
which requires only that the wave amplitude be small in 
comparison to the wavelength. All parameters are non¬ 
dimen sionalized and the analysis is made without restriction as 
to the type of velocity profile. 

Various velocity profiles are investigated and, in general, 
the phase relations between the pressure at the interface and 
wave elevation imply a situation which could allow an energy 
input from the air flow to the wave. 
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1. Intrcxiuction 


This investigation concerns a theroretical study of wind¬ 
generated water waves in which the energy transfer from a 
turbulent air flow to the underlying surface wave is examined, 

Various studies of wave generation have been made using 
models which are based on a laminar shear flow in the air boundary 
layer, and which neglect the effects of any turbulence that might 
exist in the flow. The papers of Benjamin t4 3 , and Miles C6 3 / 
£73 are of this type. However, the air flow over the ocean is 
certainly not laminar, and it seems reasonable that a study which 
takes into account the turbulence that actually exists in the 
natural environment would contribute to ourpver^dll understanding 
of the problem. It is not assumed that such a model would be the 
final answer to the problem. Indeed, it is doubtful that any single 
theory could account for all the complexities involved in the 
generation of waves on the sea. 

In the mathematical model presented here the effects of air- 
stream turbulence are incorporated by the use of mixing lengths. 

This approach may seem a bit out of style to the present-day 
hydrodynamicist, but nonetheless, it is felt that by using mixing 
lengths some useful "semi-quantititative" results might be obtained, 
at least. The general development as it is presented here is 
similar to that described by Benjamin £4 3 with the important 
exception, of course, that the viscous stresses that appear in his 
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model are replaced in this case by the Reynolds stresses associated 


with turbulent flow. 

The initial work in the formulation of this model was done by 
T. Green III (personal communication). The turbulent equations 
of motion are the results of his efforts, as are the differential 
equations which provide the basis for the numerical analysis 
described in this paper. A description of his work is therefore 
included. 

The model under consideration here is, of course, a gross 
simplification of the actual mechanism of wind-wave generation, 
and accordingly, we would expect that the results are applicable 
only under very special circumstances. It would be interesting 
to verify these results by experimentation, but such an undertaking 
is beyond the scope of this work. 

This investigation principally involves the numerical solution 
of the differential equations derived by Green, and the analysis 
of the results thereof to determine: 

1) the phase relations between the pressure at the interface 
and wave elevation, i.e., the pressure distribution over the wave 
form, and 

2) the energy transfer between the main flow and the wave-. 

2. Fonnulation of the Problem 

The air flow is assumed to be inviscid and incompressible. 

We can therefore eliminate viscous stresses as a means of energy 
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transfer, and instead we will account for the transfer of energy 
across the turbulent boundary layer by using mixing lengths. 

In mixing length theory the momentum exchange in a turbulent 
flow is related to the transverse motions of "lumps" of fluid which 
are considered to retain their identity for a certain distance L. 

This exchange of momentum can be expressed as 



where il is the mean velocity parallel to the surface, y is distance 
from the surface, the density of the air and 't' the shearing 
stress or Reynolds;Stress C 3 3 • reformulating the Navier-Stokes 
equations for turbulent flow it then seems reasonable to use Reynolds 
stresses in place of the usual viscous stress terms. Using a two- 
dimensional wave, we shall investigate only the pressure on the 
surface, and will neglect the tangential stress component. This is 
in keeping with our interest in the transfer of energy through the 
distribution of pressure over the wave, and the assumed invisdid 
properties of the fluid. 

All variables will be non-dimensionalized; the basic length 

will be the thickness of the boundary layer, D, and the basic 

velocity the speed of the flow in the region above the boundary 

layer, Thus, all lengths and velocities become multiples of 

D and Uj. Time then becomes dimensionless as D/Uj, and 

2 

stresses are expressed as multiples of 
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In establishing a frame of reference it is convenient to choose 


one in which the wave form is stationary. Thus, in a cartesian 
system of coordinates the axes will be moving with the wave at 
a speed c (the wave cfelerity), and the velocity in the flow above 
the wave in the direction parallel to the x-axis is U-c, where U 
is a function of y and represents the velocity relative to the water 


mass itself, which is moving to the left with the velocity c (not 
relative to the stationary wave form), 



Following Benjamin [4 3 / the equation of the wave is taken 
to be 

/o o\ ikx 

(2.2) yg = ae 

with the understanding that the real part represents the interface 
between the water and the atmosphere, where y is the elevation 

o 

of the surface and k is the wave number. We assume the amplitude 
a to be small in comparison with the wave length. This wave could 
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be one of a number of Fourier components which constitute the 
actual wave train. 

In order to simplify the problem and to allow for a linearized 
solution it becomes advantageous to change our frame of reference 
from the cartesian coordinates described above to an orthogonal 
curvilinear system which follows the wave form, as shown in 
Figure 2. 





Figure 2. 


The formal transformation is given by 


-k(y-ix) 


(2.3) 5 = x-iae 



Note that W = 0 now describes the interface to the order (a). 


Neglecting the normal components, the Reynolds stress in 
the new coordinate system is 


(2.5) 



(2.6) where 
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and V is a stream function to be more explicitly defined later. 

Here it should be noted that although (5 » 
intrinsic coordinates, they do approximate an intrinsic set to the 
order (^a). We assume for the purpose of calculating th6 Reynolds 
stress that the ^ coordinates are natural coordinates. Assuming 
that the fluid "lumps" of the mixing length theory move normal to 
the actual streamlines, the change in coordinate systems 

2 

introduces an error in their direction of motion of the order (a ), 

2 

and an error in the distance they travel of the order (a ). 

This implies that the stresses can be taken to be in the 
^ -direction instead of in the s-direction along the actual V-lines. 
This is shown in Figure 3. 


0 ( 3 ^) 




The mixing length relation in our curvilinear system takes 
the form L=Kd, where diis the distance between a point and the 
surface and is given by 


(2.7) 



+ a cos k^ (e"^*] -1) + O(a^) . 
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Hence, the approximation || = d is good close to the surface, but 
becomes less valid with increasing , the error bfeing of the 
order (a). Since a is assumed to be small, and the constant K is 
not firmly established, we will assume the relation L = to 
hold over the entire region 0 S ^ 1. 

The Jacobian of the transformation from cartesian coordinates 
to the curvilinear system is 



= 1 + 2kae 


“k(y-ix) 


wHlfch to the first order can be written 


(2.8) J = 1 + 2kae"’^^*] . 


Here we note that 5 ^nd are the same as the velocity 
potential and the stream function for irrotational wave motion in 
an inviscid fluid. With this in mind, let be the stream function 
in this case, its form in the absence of waves bfeing 

(2.9) * 

0 

When waves are present y may be expressed as plus a periodic 
perturbation which can be written 

(2.10) v(5,'^)“ 

where the function F has yet to be determined. 
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Following Benjamin C43 , we now formulate the appropriate 
dynamical equations in terms of ^ : 

(2.11) J(v,y„ t i J, 

(2.12) JC-'HiVy+vjru = , 

where the subscripts denote partial differentiations. 

The task now is to solve these basic equations for the 
function F and its derivatives, and finally to calculate the pressure 

I 

at the interface. 

Before proceeding further, however, it may be beneficial to 
summarize the assumptions we have made to this point. 

1) the wave amplitude is considered small. 

2) .wave length and wave number are of the order (1). 

3) the curvilinear coordinate system closely approximates a 
natural coordinate system. 

4) the mixing length L is of the form L = K f| , where K is 
taken to be 0.4. 

5) the normal component of the Reynolds stress is negligible. 

3. The Inunction F. 

We proceed in solving for the function F by first eliminating 
the pressure terms through cfsoss-differentiation and then 
combining the cross-differentiated equations to get a single 
equation in y . Linearizing by dropping terms of the order (a ) 
and higher results in 
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(3.1) ■gi[v.)V!5i *’ ■ I'jV'in ■•■ j 

+ v,v„{ 2 J, -n'JuJ 

+ + Jnj 

+ Zrj’-J^{Vw '^1)1 H'-in-i) - •^oVnoywj • 

Referring to the relations (2.8), (2.9) and 2.10), we find the 
varioCis partial derivatives of J and V with respect to 5 and 1^ , 

and substitute these into (3.1). The result is an equation in F, 

o 

which after considerable manipulation and the dropping of (a^) 
and higher-order terms can be reduced to 

(3.2) ^{(U-c)(F“-K"F)-tl"F} = rj'-ur"' 

+ Zr| r"‘(2 (JrjU" j + F'‘( 2 U'■*■ 41] 1/> r|»Lr-M|V U') 

+ e*' -4i|K +.i>|’^K’-) {u"^+ U'u")(4f| '2t|‘lc) 

+ »|VU‘''+3f U“U’"}. 

This is a fourth-order differential equation for the function F. 
To facilitate the solution we rewrite this equation as four first- 


order equations. 
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Let 


F’ = R 


F"= S 


pill— ip 


We then have 


(3.3 a) 


F' = R 


b) 


R' = S 


c) 


S' = T 


d) 


T' =F^^(R,S,T). 


Substituting in (3.2) and solving for T', i.e., F^^, gives 


(3.4) T’=i^.[K‘K{(U-c)(J-K*F)-U“F} 

-.iijT(ZU'*ilU'') -s[(Z*Ti^k‘)u‘ +4»|U" 
+ n^U”‘} - fe‘*l{u‘0”(2-4Kn +ZHV) 


This system, (3.3.a-d), lends itself quite readily to 
numerical solution. Four boundary conditions will be required to 
solve the fourth-order equation. First, the boundaries themselves 
must be established. Since the flow above the boundary layer is 
assumed to be uniform, the top of the boundary layer, i.e., |f| — 1, 
will be taken as the upper limit to the problem. This is not 
necessarily correct, but it seems as reasonable a point as any 
other to use, and computations are simplified if we do not encounter 
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the transition from a boundary-layer velocity profile to a uniform 


velocity in the region above = 1. 

Since our model is concerned only with turbulent flow, the 
laminar sublayer is neglected. Accordingly, we will take the top 
of the laminar sublayer to be the lower boundary, assuming its 
thickness to be on the order of 0.01 to 0.1 cm. 

The velocity components of the flow parallel to the ^ ^nd 
axes respectively are 


(3.6) n= . U-c ♦ a(F'+uV''')e'” 

(3.7) v= J‘‘v, = -lKa{F- + (W-c)e'"'}a“‘J. 

X 

At the lower boundary, the velocity U is zero. The normal 
component v must also be zero since the wave is stationary; 
hence (3.7) gives 


(3.8) F(0) = c. 

The velocities in the water must take on negative values in order 

that the shearing stress due to the air flow can be continuous 

across the interface. We assume that the variable component of 

the tangential velocity is expressed in the form Bae^^S , where 

B = —. The tangential velocity then becomes u = -c-Bae^^^ , 

c 

and from (3.6); 

(3.9) ‘ F' = -U'(0) + B. 
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Turning our attention to the upper boundary, we again find 
the normal velocity v to be zero since the flow is assumed 
parallel to the f|-lines. Here U = 1, and (3.7) gives 

(3.10) F = (c-l)e“’^*I . 

U is assumed to be a constant with the value 1.0 at the upper 
boundary and above, hence U' becomes zero and the tangential 
velocity u = 1-c. Applying these conditions to (3.6) yields 

(3.11) F' = 0. 

These boundary conditions are consistent with the linearized 
expression for y (2.10), and their validity depends only on ka 
being small. 

We now proceed with the solution for the function F. Several 
numerical methods are available for solving the system of linear 
differential equations (3.3a,b,c) and (3.4). However, since 
each method involves an iterative procedure that begins at one 
boundary and ends at the other, it is not possible to include all 
the boundary conditions as initial values. Therefore, a single, 
direct solution to the set of equations is not feasible. Instead, 
five separate solutions will be found, each starting at the same 
boundary with slightly varied initial conditions, and the final 
solution will be obtained by superposing the results to fit the 
required conditions at the other boundary. 

Using the upper boundary as a starting point, we pick the 
initial conditions for the first four solutions in the form 


14 
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(3.12a) 


F^(l) = (1) = f;- (1) = 0 and F^" (1) = 1 

b) F^d) = F-^d) = F^'d) = 0 and F^’d) = i 

6) F^d) =F^d) =Fyd) = 0 and F^'(l) = 1 

d) F^d) =F^d) =F-'(1) =OandF-(l) =i . 

Using these initial conditions and the homogeneous form 
of (3.4), which is 


( 3 . 13 ) T‘» ipu.[K*K{(o-c)(^-K‘r)-u"r} 
-ZijTfZUVijy") - ^ {(if r|*-K*)U'+Aij U "f ,|»U"*}] 


we evaluate F and F' at the lower boundary, 
have the form 


(3.14a) 

Fj^(O) - a^^ + ib^ 

F'j(O) =c 

b) 

*'=2 

F^(0) = c 

c> 

FjlO) = aj + ibj 

F^(0) = c 

d) 

^ 4 ( 0 ) = ^4 + ‘b4 

F'(0) =c 
4 


The four results 


+ id. 


+ id. 


+ id. 


+ id.. 
4 


For the fifth solution we use the prescribed upper boundary 
conditions as initial conditions 

(3.12e) . Fgd) = e’’^(c-l) and F^(l) = F” (1) = F^"(l) = 0. 
and again evaluate at the lower boundary, this time using the 
inhomogeneous equation (3.4). The result is 
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F^(0) = O5+ld5. 


(3.14e) F5(0)=a5 + ib5 


We can now superpose the five values found for F and F' 
at the lower boundary to match the prescribed lower boundary 
conditions. Let 


(3.15) 

^ = Vl V2 * ^ 3’’3 V4 '■5 

(3.16) 



where the J's are real constants. The lower boundary conditions 
now prescribe the equations 


(3.17^) 

Ta +Ja„+Ta +J,a +a =c 

11 2 2 ^3 3 4 4 5 

b) 

^ ^ ^ '=5 = “ 

c) 


d) 

" Vs " '4^4 ^ ^ ° 


which form a set of simultaneous linear equations that can be 
readily solved for the constants through , 

Using the expressions (3.15) and (3.16) we can find the 
values of F and F' at the lower boundary. The values of F” 
and F'" can be found Similarly. 

3. Pressure 

We assume that the energy transfer between the air and the 
surface wave is accomplished through the normal pressure 
component at the surface. An estimate of the pressure 
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distribution over the wave is therefore necessary. As in 
Benjamin C43 , we consider the pressure variation due to the 
wave disturbance as distinct from that in the primary flow. We 
then define this "perturbation pressure" as 

(4.1) p = P(H)ae^’'^ 

where P is a function of only. Differentiating with respect 
to ^ : 


(4.2) p^ = ikp. 

Using the relations (4.1) and (4.2) with 2.11) and solving 
for p gives 

(4.3) P(H ) = U'Fa(U-c)F' - 2iK2 fj (^ *J2U'(F" + U"e"’ 

+ rje“^1u" [u'-k2(u-c)j + rjU'(F"'U'"e" 

Alternatively, we can solve for P by using (2,12) and 
integrating p^^ with respect to r| , making use of the fact that 
the disturbance of the flow vanishes at some distance from the 
surface. To be precise, the vanishing point should be taken 
as infinity. However, for this study we will assume that the 
perturbation becomes negligibly small at the upper limit of the 
boundary layer, i.e., at = 1, and this will be used as the 
upper limit for integration. P in this case has the form 


M.4) 




.)Fa» 


2.iKV[ + +r|kU*F|an . 
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The pressure at the surface could be obtained from either 


of these expressions, which obviously must be equivalent. For 
ease of computation only (4.3) will be used. 

The perturbation pressure solution will be in the form 
P( r| ) = P^+ P^, where the imaginary part P^ represents the 
component out of phase with the wave form, with a negative 
value indicating a phase shift to the left in the case of a wave 
moving to the right. This would produce a lagging pressure 
perturbation and provide a situation in which energy is 
transferred to the wave (see below). 

5. Energy Relations 

We now consider the energy transferred from the air flow to 
the wave, the energy dissipated within the wave, and the effect 
of the net energy input into the wave in terms of wave growth. 

The energy transfer between the two media, as has been 
pointed out, is assumed to be dependent solely upon the 
distribution ,6f pressure over the wave. Specifically, it is 
dependent upon the phase relations between the perturbation 
pressure and the wave form. If they are in phase there is no 
energy transfer. Maximum interchange of energy will occur 
when the phase difference is 90 degrees, with an energy input 
to the wave taking place when the pressure lags behind the 
wave. This relation places the area of highest pressure on the 
upwind side of the wave crest, and looking at the mechanics of 
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the situation intuitively, we can think of the wind as pushing 
the wave from behind. 

The rate of energy transfer from the air flow to the wave per 
unit surface area is 



On integrating and dimensionalizing, this becomes 

(5.2) Eg= la2kP^(j)^U^). 

The rate of energy dissipation due to the viscosity of the 
water is 

(5.3) EL=2jUa2k^c2 

where ^ is the coefficient of viscosity of the water, 
with a value of approximately 0.015 dyne-sec/cm^. 

Subtracting (5.2) from (5.1) we are able to find the net 
energy added to the wave per unit surface area per unit time, 
which is the energy available for wave growth. Using the 
relation 

and differentiating with respect H gives an expression for the 
rate of wave growth in the form 

4 

(5.4) dH= -—En. 
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The wave growth per cycle (i. e., the amount the wave will 
grow during the time required to travel one wave length) is 
given by 


(5.5) 


dH = 


8 E 


n 




where the term D/U^ is used to dimensionalize the non- 
dimensional c and k. 

6. The Computer Solution 

The problem was solved numerically on a CDC 1604 digital 
computer using Fortran 63 as the programming language. The 
computer program is included a s Appendix I. 

The program is designed to obtain a solution following the 
procedures outlined in Sections 3,4, and 5. Pressure at the 
Interface is computed, as i^energy exchange in the wave. 
Percentage wave growth per wave length is also evaluated in 
order to give a more meaningful end result. 

The system of linear differential equations (3.3) was solved 
using a subroutine which employs the Runge-Kutta Gill Fourth- 
Order method C 2 J . This method is an iterative procedure, 
the accuracy of which depends on the size of the increment of 
that is used. The process was started in each case at the upper 
boundary (|^ = 1). An initial increment size of 0.001 was found 
to be the optimum. An increase to 0.01 produced results which 
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were considerably different, whereas reducing the size to 0.00025 
gave only a two percent "improvement" in the evaluation. In the 
region below = 0.001 the increment size was reduced as 
necessary to reach the lower boundary point. 

The set of simultaneous linear equations (3.17) was solved 
with a second subroutine which employs a modification of the 
Jordan Elimination method jQlJ . Both subroutines were 
obtained from the files of the computer facility at the U. S. Naval 
Postgraduate School. 

7. Results and Conclusions 

The results of this study are presented as evidence of the 
applicability of the method of solution and not necessarily as 
proof of the validity of the model. 

It was initially intended to examine a series of velocity 
profiles at various wave speeds and to investigate the Importance 
of the assumed thickness of the laminar sublayer by varying the 
location of the lower boundary. However, due to the unexpectedly 
long solution time (12 to 13 minutes of computer operation) the 
investigation was limited to four velocity profiles, only one of 
which was evaluated at more than one lower boundary location. 

The four velocity profiles used were df the forms: 

(7.1) U = 1 + w ln( ) where w = -1/In(lower boundary) 

(7.2) t ; U = H 
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(7.3) 


U=sin(l/ltfr) ) 


(7.4) U = l-(1-K|)^. 

An arbitrary boundary layer height (D) of ten meters was 
used, and the thickness of the laminar sublayer was t^en to 
be 0.01 cm. This placed the non-dimensionalized lower 
boundary at = 1x10 ^. The wind speed above the boundary 
layer (Uj^) was chosen as ten meters per second. Each profile 
was treated separately using values of the non-dimensional wave 
speed c ranging from 0.0 to 1.0. 

The logarithmic profile (7.1) and the exponential profile 
(7.2) were only partially examined. Both were abandoned after 
the first few runs when it appeared that neither would produce 
meaningful results. In both cases the values for Surface 
pressure and wave growth rate were of unreasonable magnitudes. 
This is attributed to the extreme values of the derivatives of 
these particular profiles in the region close to the lower 
boundary. Table 1 summarizes the results obtained using these 
two profiles. 

The quarter-period sinusoid profile (7.3) and the cubic 
profile (7.4) produced results that seemed reasonable and both 
were examined over the entire c-range (0.0 to 1.0). These 
results are shown in Tables 2 and 3, and Figures 4 and 5. 
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In the case of the sinusoid profile, a damping of the wave 
was indicated at the higher wave speeds. Wave growth was 
indicated between about c = 0.35 and 0.19, with a maximum 
growth rate occurring near c = 0.2. Damping was again evident 
at c = 0.18, but below this point the growth rate increased 
sharply. The region between 0.15 and 0.1 was not examined 
and in view of the large discrepancy between the values obtained 
at these two points (a large growth rate at 0.15 and extreme 
damping at 0.1) there is apparently no predictable trend in this 
region, 

The cubic profile indicated wave damping throughout the 
entire c-range, with maximum damping occurring between c = 0.5 
and 0.6. As c was reduced below 0.2 the pressure oscillated 
with increasing amplitude and frequency, producing alternating 
values of extreme wave growth and severe wave damping. 

The extreme oscillations in wave growth at the lower wave 
speeds would tend to indicate that the model is not applicable 
in this region. From the data obtained it appears that the 
behavior of the model in this range of wave speeds is affected 
by both the velocity profile and the thickness of the laminar 
sublayer. High growth rates, in themselves, are not to be 
entirely discounted since such phenomena have been observed 
in the Initial stages of wave development. This nearly 
instantaneous wave growth is mentioned by Gelci et al. C 5 3 
and Walden C 8 3 and C 9 3 • 
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The quarter-period sinusoid profile was subsequently 
examined using 0.05 cm, as the thickness of the laminar 
sublayer. This placed the lower boundary at = SkIO"^. It 
is interesting that the results were considerably different from 
those obtained with = 1x10 ^ as the lower boundary. Wave 
damping was indicated over the entire range and varied inversely 
with the wave speed, the damping rate increasing sharply as c 
was decreased below 0.2. The results are shown in Table 4 and 
Figures 4 and 5 . 

Further investigation of the model can be readily accomplished: 
additional wind velocity profile can be used, the boundary layer 
height (D) and the wind speed above the boundary layer (U^^) can 
be varied, as can the thickness of the laminar sublayer. In 
addition, the model itself could be altered by a re-derivation 
used a different form of the mixing length. Perhaps a mixing 
length based on actual distance from the surface (2,3) rather 
than QnL= kt| may prove to be a more valid form. 

Although the data obtained was not as extensive as had 
been planned, the following conclusions can be drawn 
concerning the model: 

1) The velocity profile is critically important; those profiles 
having derivatives of high orders of magnitude near f| = 0.0 
cannot be used. 
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2) The results are highly dependent upon the location of 
the lower boundary. 

3) The model may not apply at the lower wave speeds; its 
applicability in this speed range apparently depends upon the 
velocity profile and the laminar sublayer thickness. 

4) Further investigation of the model is warranted. In 
particular, more numerical solutions of this (dL = k|^) model should 
be obtained. Then a more reasonable form of the mixing length 
should be Introduced (see above), although this will complicate 
the resulting differential equation for F. 

In conclusion, the results^bbthlned are very different from 
those obtained by Miles C63 . Although this obviously does 
not negate that work, it does suggest that the effect of turbulence 
should not be heglected in the study of wind-generated water 

I 

waves. 
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c 

Wave 

Length 

(Meters) 

Wave Growth 
Per Cycle 
(Percent) 

Perturbation 
at the Laminar 
Real 

Pressure 

Sublayer 

Imag 

.100 

.64 

-7975009.35 

5409.0000 

53160.0000 

.200 

2.56 

-18939896.02 

7780.0000 

252500.0000 

.300 

5.77 

-27213442.61 

9155.0000 

544200.0000 

.500 

16.02 

2219980.83 

-5369.0000 

-73990.0000 



Logarithmic Velocity Profile 


. 100 

.64 

-494.94 

-4.3220 

3.2960 

. 150 

1.44 

-948.66 

-11.3400 

9.4840 

.300 

5.77 

-2942.89 

-44.2400 

58.8500 

.500 

16.02 

-6036.77 

-103.9000 

201.2000 



Exponential Velocity Profile 



Table 1. Non-dimensional perturbation pressure and percentage 
wave growth per cycle as a function of non-dimensional 
wave cererity, c, and wave length for logarithmic and 
exponential velocity profile. (The thickness of the 
laminar sublayer is taken as 0.01 cm, the boundary 
layer-height as 10 meters and the wind speed above 
the boundary layer as 10 meters per second.) 
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c 

Wave 

Length 

(Meters) 

Wave Growth 
Per Cycle 
(Percent) 

Perturbation 
at the Laminar 
Real 

Pressure 

Sublayer 

Imag 

.100 

.64 

-41495.73 

140.8000 

276.6000 

.150 

1.44 

293.99 

-.2273 

-2.9410 

.160 

1.64 

67,69 

-1.1030 

-.7232 

.170 

1.85 

14.85 

-.5639 

-.1694 

.180 

2.08 

-3.08 

-.6428 

.0360 

,190 

2.31 

.18 

-.7120 

-.0032 

.200 

2.56 

3,32 

-.6778 

-.0450 

.300 

5.77 

1.59 

-.5443 

-.0321 

.350 

7.85 

.45 

-.4801 

-.0108 

.400 

10.25 

-.58 

-.4097 

.0153 

.500 

16.02 

-4.37 

-.2837 

.1456 

.600 

23.07 

-3.25 

-.2593 

. 1298 

.700 

31.39 

-2.81 

-.1828 

.1309 

.800 

41.01 

-2.54 

-.1537 

.1355 

.900 

51.90 

-2.33 

-.1562 

. 1400 

1.000 

64.07 

-2.20 

-.1778 

.1469 


Table 2. Non-dimensional perturbation pressure and percentage 
wave growth per cycle as a function of non-dimensional 
wave cererity, c«and wave length for a quarter-period 
sinusoid velocity profile. (The thickness of the 
laminar sublayer is taken as 0.01 cm, the boundary 
layer height as 10 meters and the wind speed above 
the boundary layer as 10 meters per second.) 
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c 

Wave 

Length 

(Meters) 

Wave Growth 
Per Cycle 
(Percent) 

Perturbation 
at the Laminar 
Real 

Pressure 

Sublayer 

Imag 

.100 

.64 

22217.33 

-11060.0000 

-148.1000 

, 150 

1,44 

-2678.48 

44.7400 

26.7800 

.160 

1,64 

-354.63 

-1.5610 

3.7810 

. 165 

1.74 

-379.25 

3,9680 

4.1700 

.170 

1.85 

135.01 

-.5818 

-1.5310 

. 180 

2.08 

77.21 

-.8456 

-.9274 

.190 

2,31 

-395,17 

1.7300 

5.0040 

,200 

2.56 

-3.39 

-.6039 

.0444 

.300 

5,77 

-5.26 

-.2775 

. 1048 

.400 

10.25 

-7.89 

-.0582 

,2102 

.500 

16.02 

-11.61 

.0199 

.3869 

.600 

23.07 

-12.47 

-.0833 

.4987 

. 700 

31.39 

-9.14 

-.1308 

.4263 

.800 

41.01 

-7.40 

-.0985 

.3947 

.900 

51.90 

-6,75 

-.0757 

.4049 

1.000 

64.07 

-6.51 

-.0696 

.4341 


Table 3. Non-dimensional perturbation pressure and percentage 
wave growth per cycle as a function of non-dimensional 
wave cererity, c^and wave length for a cubic velocity 
profile, (The thickness of the laminar sublayer is taken 
as 0.01 cm, the boundary layer height as 10 meters and 
the wind speed above the boundary layer as 10 meters 
per second.) 
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c 

Wave 

Length 

(Meters) 

Wave Growth 
Per Cycle 
(Percent) 

Perturbation 
at the Laminar 
Real 

Pressure 

Sublayer 

Imag 

o 

o 

1—1 

.64 

-117.65 

1.0800 

.7810 

.150 

1.44 

-23.56 

.3071 

.2341 

.160 

1.64 

-20.65 

.3315 

.2190 

.170 

1.85 

-19.48 

.3367 

.2196 

.180 

2.08 

-18.66 

.3469 

.2229 

.190 

2.31 

-17.75 

.3554 

.2239 

.200 

2.56 

-16.88 

.3620 

.2242 

.300 

5.77 

-12.34 

.4338 

.2464 

.400 

10.25 

-10.21 

.5092 

.2721 

.500 

16.02 

-9.67 

.5941 

.3222 

.600 

23.07 

-8.56 

.6136 

.3424 

.700 

31.39 

-7.56 

.6517 

.3526 

,800 

41.01 

-6.98 

.6746 

.3723 

.900 

51,90 

-6.69 

.6871 

.4012 

1/000 

64.07 

-6.53 

.6942 

.4352 


Table 4. Non-dimensional perturbation pressure and percentage 
wave growth per cycle as a function of non-dimensional 
wave cererity, c, and wave length for a quarter-period 
sinusoid velocity profile. (The thickness of the laminar 
'%!W«fiublayer is taken as 0.05 cm, the boundary layer height 
as 10 meters and the wind speed above the boundary 
layer as 10 meters per second.) 
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Figure 4. Percentage wave growth per cycle vs. non-dimensional wave celerity c. 













Figures. Out-of-phase non-dimensional pressure amplitude P. vs. non-dimensional wave celerity, c. 
A negative value of Pj indicates energy transfer to the wave. 
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' appendix I 

* 0 
rl ^ 

THE COMPUTER PROGRAM. 

THIS PROGRAM EMPLOYS FORTRAN 63 LANGUAGE AND IS DESIGNED FOR 
USE WITH THE CDC 1604 COMPUTER, 


PROGRAM WIMBO 

TYPE COMPLEX Y,F,YF,YIN,FCNF♦VKAPA,VKAPA2»P,PN»PG♦BIGP.H.WK,U,U1» 
1U2»U3'.U4 

DIMENSION Y(4»5)»FCNG(2»4)»G(4»5),CONST(4)»FCNF(4)»YF(4,5)»F(4) 
w 1,YIN(4) 

C 

C WAVE CELERITY INPUT 

C=.8 

' C 

DO 71 N2=l,5 
DO 71 Nl=1.4 
71 Y{N1»N2)={0.»0.) 

Y<4,1)={1.»0.) 

, V Y { 4,2) = (0 . , 1. ) 

Y(3»'3‘r3'(l.»0.) 

Y(3,4)=(0.,l.) 

C 

C BOUNDARY LAYER HEIGHT INPUT (CM) 

DELTA=1000. 

C 

C WIND SPEED INPUT (CM/SEC) 

VACT=1000. 

C 


C 


BETA=980,665*DELTA/(C»VACT**2) 







nn non nnnn 


) 


C WAVE NUMBER (NON-DIMENSIONALIZED) 

WK=980.665*DELTA/(C*VACT )**2 


WAVE LENGTH (METERS) 

WL=2.*3.14159/WK*DELTA/100, 

WV=VACT/100. 

BL=DELTA/100. 

PRINT 390 
390 FORMAT (IHl) 

PRINT 39*C»WK»BETA*8L*WV»WL 

39 FORMAT (28HOSINUSOIDAL VELOCITY PROFILE, 5X,3HC= F4.2,5X,3HK= F6.2 
1,5X»6HBETA= F5.4,5X,25HBOUNDARY LAYER THICKNESS= F3.0»7H METERS/ 
212H0WIND SPEED= F3.0»11H METERS/SEC 5X,12HWAVE LENGTH= F6.1,7H MET 
. 3ERS) 

Y(1,5)=EXPF(-WK)*(C-1.) 

DO 88 J=l,5 
ETA=1. 

DO 66 IN=1,4 
66 YIN(IN)=Y(IN»J) 

PRINT 3, J,(YIN(N)♦N=l,4) 

3 FORMAT (/////30H0BOUNDARY CONDITIONS INPUT J = 12,15X,4C(FIO.1»F5. 
11 ) ) 

EVALUATE THE FUNCTION F 

CALL FUNCTF (YIN,J,WK,C»FCNG,U1»ETA,U2,U3,U) 

DO 75 13=1,4 
75 YF(13,J)=YIN(13) 

IF (J-5) 21,22,21 
21 CONTINUE 

matrix INPUT FOR THE FIRST FOUR SOLUTIONS (HOMOGENEOUS) 

G(1,J)=FCNG(1,1) 

G(2,J)=FCNG(2,1) 
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G(3»J)=FCNG{1,2) 

G(4,J)=FCNG(2»2) 

GO TO 23 

MATRIX INPUT FOR THE LAST SOLUTION (INHOMOGENEOUS) 

.22 G(1»J)=-FCNG(1,1)+C 
G(2»J)=-FCNG(2»1) 

G(3»J)=-FCNG(1,2)~U1+BETA 
G(4,J)=-FCNG(2»2) 

23 CONTINUE 

PRINT 31* (G(I,J)♦1=1,4)» (YIN(L),L=1,4) 

31 FORMAT (13HOMATRIX INPUT//15X,4E12.3*/20H0SUBFUNCTION F INPUT//15X 
1♦4C{E12.3,E12.3)) 

88 CONTINUE 
PRINT 34 

34 FORMAT (IIHOTHE MATRIX/) 

PRINT 35* ( (Gdl.JJ) *JJ=1*5) *I 1 = 1*4) 

35 FORMAT {5E15.5/) 


\ 


j 

{ 

I 

i 

I 

♦ 

1 


.A 


SOLVE FOR THE CONSTANTS J1 THROUGH J4 
CALL WIM80Z0 (G*CONST) 

PRINT 38 

38 FORMAT (/////28HOTHE FUNCTIONS FI THROUGH F5 79X*13HTHE CONSTANTS 
- r/107X*13HJl THROUGH J4 //) 

PRINT 32* ((YF(I5*J1)* 15=1*4)*CONST(J1)*Jl = l*4)*(YF(16*5)*16 = 1*4) 
32 FORMAT (/4(4C(E14.3*E12.3),E15,3//)*4C{E14.3 *E12.3)) 

DO 76 14=1,4 
F(I4)=0. 

DO 77 M=l*4 

77 F(I4)=F(I4)+C0NST(M)*YF(I4*M) 

F(I4)=F(I4)+YF{14*5) . • . ^ 

76 CONTINUE 
PRINT 390 

PRINT 39* C*WK*BETA*BL*WV*WL 
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PRINT 33» ETA»(F(M2)»M2=1»4) 

33 FORMAT (//30H0VALUE OF FUNCTIONS F AT ETA = fFS.6//5X»4C(El5.3,E12 
1.3) ) 

ETAK=ETA*Wk 

ENK=EXPF(-ETAK) 

~ WK2=WK*WK 

VKAPA2=(0.,0.32) 

SURFACE PRESSURE 

BIGP=U1*F(1)-(U-C)*F(2)-VKAPA2*ETA/WK*(2.*U1*(F(3)+U2*ENK)+ETA»ENK 
1»U2»(U1-WK2*(U-C))+ETA»U1*(F(4)+U3*ENK)) 

PRINT 40» BIGP 

40 FORMAT {//30H0THE VALUE P COMPUTED FROM-PXI »C(E15.3,E15.3)) 

CALCULATE ENERGY RELATIONS 

CALL ENERGY (C,WK»BIGP,DELTA,VACT) 

100 C=C+.3 
99 END 


SUBROUTINE FUNCTF (Y»JtWK»C,FCNF»U1,ETA»U2»U3jU) 

THIS SUBROUTINE EVALUATES THE FUNCTION F AND ITS DERIVATIVES AT THE 
TOP OF THE laminar SUBLAYER 

" TYPE COMPLEX Y,F,VKAPA,VKAPA2,FCNF,P,PN»PG,BI6P,H,K,U,01,U2♦U3,U4, 
1W<,WK2»ETAK,ETAK2»ETA2U»ETA3»ENK»ETA2 
DIMENSION Y(4)» F(4)» FCNF(4) 

DATA (NT=0), (N=4)» (VKAPA=(0,,3.125)) 

PRINT 39» (YIN),N=1»4) 

39 FORMAT (47X,4C(FIO.1,F5.1)) 

PI=3.14159/2. 

CL=-l./LOGF(.00001) 

11 CONTINUE 

IF (ETA-.0001) 80»80»81 
80 H=-.00001 
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GO TO 200 

81 IF (ETA-.OOl) 82*82,83 

82 H=-.0001 
GO TO 200 

83 H=-,001 
200 CONTINUE 

EP=PI*ETA $ U=SINF{EP) S U1=PI*COSF(EP) $ U2=-PI*PI*SINF(EP) 
U3=-(PI**3)*COSF(EP) $ U4=PI**4*SINF(EP) 

ETA2=ETA**2 
ETAK=ETA*WK 
ETAK2=ETAK**2 
.WK2=WK*WK 
ETA2U=1/ETA2/U1 
F(1)=Y(2) 

F(2)=Y(3) 

— F(3)=Y(4) 

IF (J-5) 600,601,600 
C 

C HOMOGENEOUS EQUATION 

600 F(4)=ETA2U*(WK*VKAPA*(tU-C)*(Y(3)-WK2*Y(1))-U2*Y(1))-2.*ETA*Y(4)*( 
12.*U1+ETA*U2)-Y(3)*((2.+ETAK2)*U1+4.*ETA*U2+ETA2*U3)) 

. GO TO 603 
C 
C 

C INHOMOGENEOUS EQUATION 

. 60l"Fr^T*=m2U*(WK*VKAPA*( (U-C )»(Y ( 3 )-WK2*Y (1) )-U2*Y( 1) )-2.*ETA*Y( 4 )* ( 
12.*U1+ETA*U2)-Y(3)*((2•+ETAK2)*U1+4.*ETA*U2+ETA2*U3)-EXPF(-ETAK)*( 
2U1*U2*(2.-4.«ETAK+2.*ETAK2)+(U2**2+U1*U3)*(4.*ETA-2.*ETA2*WK)+ETA2 
3*U1*U4+3.*ETA2*U2*U3)) 

C 

603 CONTINUE 
C 

C SOLVE THE SYSTEM OF LINEAR DIFFERENTIAL EQUATIONS 

S=CRKLDEQ(N,Y,F,ET A,H,NT) 

IF(S-1.0)99,11,12 
12 CONTINUE 
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C LOWER BOUNDARY CUTOFF 

IF (ETA-.000016) 99,99»11 
C 

99 CONTINUE 
DO 61 I=l»4 
61 FCNF(I)=Y(I) 

PRINT 30, ETA, U, (FCNF(I),I=1,4) 

30 FORMAT (2IHOSUBFUNCTI ON F OUTPUT //F9.6,4X,E9.2,4C(E12.3,E12.3) ) 
RETURN 
END 


FUNCTION CRKLDEQ (N,Y,F,X,H,NT) 00000 

THIS-SUBROUTINE SOLVES THE SYSTEM OF LINEAR DIFFERENTIAL EQUATIONS 
BY THE RUNGE-KUTTA GILL METHOD 

DIMENSION Y(10),F(10),Q(10) 00010 

TYPE COMPLEX Y,F,Q,X,H ‘ 00020 

NT=NT+1 00030 

GO TO (1,2,3,4),NT 00040 

I DO 11 J=1,N 00050 

II 0(J)=0, 00060 

A=.5 00070 

X=X+H/2. . 00080 

GO TO 5 00090 

2.- AP=.7g'289321881 00100 

GO TO 5 00110 

3 A=l.7071067812 00120 

X=X+H/2. 00130 

GO TO 5 00140 

4 DO 41 1=1,N 00150 

41 Y(I)=Y( I )+H*F(I)/6.-Q( I)/3. 00160 

NT=0 00170 

CRKLDEQ=2. 00180 

GO TO 6 00190 

5 DO 51 L=1,N 00200 
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Y(L)=Y(L)+A*(H*F(L)-Q(L)) 00210 

51 Q(L)=2.»A^^H^^F(L) + ( l.-3.*A)*Q(L) 00220 

CRKLDEQ=1. 00230 

6 CONTINUE 00240 

RETURN 00250 

■ END 00260 


SUBROUTINE WIMBOZO (A»X) 

THIS SUBROUTINE SOLVES FOR THE CONSTANTS J1 THROUGH J4 USING A 
MODIFICATION OF THE JORDAN ELIMINATION METHOD 


DIMENSION A(4»5),X(4) 

K=5 $ N=4 

11 IF (K-l)13,6»15 
15 0=0* 

DO 2 I=2»K 

IF(ABSF(A(I-l,l))-D) 2»2»3 

3 L =1-1 

D = ABSF(A(L»1)) 

2 CONTINUE 

4 IF(L-1)5»6»5 

5 DO 7 J=1»K 
D =A(L»J) 

A(L»J) = A(1»J) 

- 7'Atl»J) = D 

6 DO 8 I =1»N 

8 X(I) = A( 1,1) 

IF (<-1)12,13,12 

12 DO 10 J=2,K 

D = A(1,J)/X(1) 

DO 9 1=2,N 

9 A{I-1,.J-1) = A(I,J) -X(I)*D 

10 A(N,J-1) =D , '1 

K = K~1 
GO TO 11 


00040 

00050 

00060 

00070 

00080 

00090 

00100 

00110 

00120 

00130 

00140 

00150 

00160 

00170 

00180 

00190 

00200 

00210 

00220 

00230 

00240 

00250 
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13 CONTINUE 
RETURN 
END 


00260 

00270 

00280 


i 
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ENERGY RELATIONS 

THIS SUBROUTINE CALCULATES THE VAR 


lOUS ENERGY RELATIONS 


SUBROUTINE ENERGY (C»WK»P»DELTA♦VACT) 
TYPE COMPLEX WK 
DIMENSION P(2) 


WAVE AMPLITUDE {NON-DIMENSlONALIZED) AND VISCOSITY 
A=.01 $ VISC=*01 

DENSITIES OF AIR AND SEA WATER 

RHO=.0012 

RHOW=1.025 


ENERGY TRANSFER TO THE WAVE 
EGAIN=.5*A*A«WK*{-P{2))*RHO*VACT**3 


C 


C 


ENERGY LOSS IN THE WAVE 

ELOSS=2.*VISC*A*A*WK**3/DELTA*C*C*VACT»*2 


2 


1 

3 


enej=egain-eloss 

DE=2>3.14159*DELTA*ENET/(WK*C*VACT) 

WAVE GROWTH (PERCENT) 

DA=4.*DE/(RHOW*980.665*A*DELTA)*100, 

PRINT 2 ♦EGAINjELOSS 

FORMAT (/22H0ENERGY GAIN FROM WIND E17.5/26H ENERGY DISSIPATED IN 
IWAVE E13.5) 

FORMAT^(//46H0NET ENERGY TRANSFER TO THE WAVE PER UNIT TIME E15.5) • 
PRINT 3» DA 

FORMAT (33H0PERCENT GROWTH OF WAVE PER CYCLE F7*2) 


END 


I 
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